function depvars = Projection(cfg,exvars,depvars,indepvars,r)

%% Reconstruction of historic set
modes  = indepvars.UseModes(r,1,:);
for t = 1:35
    eval([ '[junk depvars] = ' cfg.Model '(cfg,exvars,depvars,indepvars,r,t);']);
    GeneralPlotting('GUIPlots','',cfg,exvars,depvars,indepvars,r,t);
    if cfg.ParallelComputing==0; set(cfg.h.Info,'String',{'Region / Time / Run: ',[num2str(r) ' / ' num2str(t) ' / ' num2str(cfg.Run)]}); end
end
    
%% Determine errors (based on NRSME)
depvars.errors.TTBfit(r)   = sqrt(sum((sum(depvars.TravelDemand(r,:,modes)./exvars.Speed(r,:,modes),3) - ...
    indepvars.TTB(r,1)).^2./indepvars.TTB(r,1).^2,2)/size(depvars.TravelDemand(r,:,modes),3));

depvars.TMB(r,:)              = (exvars.TMBCorrection(r,1) + 0.035 + 0.085 .* sum(depvars.ModeSplit(r,:,cfg.ModesFasterThanLDV),3));
depvars.errors.TMBfit(r,1)   = sqrt(sum((sum(depvars.TravelDemand(r,:,modes).*exvars.TravelPrice(r,:,modes),3)./exvars.Income(r,:) - ...
    depvars.TMB(r,:)).^2./depvars.TMB(r,:).^2,2)/size(depvars.TravelDemand(r,:,modes),3));

depvars.errors.TDfit(r,1)    = sqrt(sum(((sum(depvars.TravelDemand(r,cfg.CalibrationTime,modes),3) - sum(exvars.TravelDemand(r,cfg.CalibrationTime,modes),3)) ...
    ./sum(exvars.TravelDemand(r,cfg.CalibrationTime,modes),3)).^2,2)/length(cfg.CalibrationTime));

depvars.errors.ModeFit(r,1)  = sqrt(sum(sum((exvars.TravelDemand(r,cfg.CalibrationTime,modes)...
    ./repmat(sum(exvars.TravelDemand(r,cfg.CalibrationTime,modes),3),[1 1 length(cfg.modes(modes))])-depvars.ModeSplit(r,cfg.CalibrationTime,modes)).^2,3))/length(cfg.CalibrationTime)); 

depvars.errors.OFfit(r,1)    = sqrt(sum(depvars.ObjectiveFunction(r,cfg.CalibrationTime).^2)/length(cfg.CalibrationTime)); 

clear junk
%% Plot the minimization output
GeneralPlotting('Visualize parameter space','',cfg,exvars,depvars,indepvars,r, cfg.CalibrationYears);
GeneralPlotting('GUIPlots','',cfg,exvars,depvars,indepvars,r, cfg.CalibrationYears);
GeneralPlotting('PlotMinimize','',cfg,exvars,depvars,indepvars,r, cfg.CalibrationYears);

%% Forecast
if strcmp(cfg.Model,'ModelPoles') || strcmp(cfg.Model,'ModelSimple')
     depvars.TravelDemand(r,35,modes) = exvars.TravelDemand(r,35,modes);
     depvars.ModeSplit(r,35,modes)    = exvars.ModeSplit(r,35,modes);
    eval([ '[junk depvars indepvars] = ' cfg.Model '(cfg,exvars,depvars,indepvars,r,36);']);
    
elseif strcmp(cfg.Model,'ModelGcam')
    eval([ '[junk depvars] = ' cfg.Model '(cfg,exvars,depvars,indepvars,r,35);']); 
    indepvars.Sigma(r)  = ...
        sum(exvars.TravelDemand(r,cfg.CalibrationYears(end),modes),3)...
        ./((exvars.Income(r,cfg.CalibrationYears(end))^indepvars.IncomeElasticity(r)).* ...
        (sum((indepvars.TravelPrice (r,cfg.CalibrationYears(end),modes)...
        + exvars.WageRate (r,cfg.CalibrationYears(end))...
        ./ exvars.Speed(r,cfg.CalibrationYears(end),modes))...
        .*depvars.ModeSplit(r,cfg.CalibrationYears(end),modes),3)^-indepvars.PriceElasticity(r)));
end

for t = 36: length(cfg.Time)
    eval([ '[junk depvars] = ' cfg.Model '(cfg,exvars,depvars,indepvars,r,t);']);
    GeneralPlotting('GUIPlots','',cfg,exvars,depvars,indepvars,r,t);
    if cfg.ParallelComputing==0; set(cfg.h.Info,'String',{'Region / Time / Run: ',[num2str(r) ' / ' num2str(t) ' / ' num2str(cfg.Run)]}); end
end


end